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1.  Introduction 


Light  detection  and  ranging  (LiDAR,  also  known  as  LADAR)  technology  emerged  soon  after  the 
introduction  of  the  laser  in  the  late  1950s,  but  it  is  only  widely  adopted  recently  for  a  myriad  of 
applications  in  various  engineering  fields,  as  well  as  for  use  in  the  geospatial  intelligence 
community  and  the  military  (/).  Many  types  of  LiDAR  systems  have  been  developed  to  date, 
such  as  two-dimensional  (2-D)  gated  framing  LIDAR,  scanning  three-dimensional  (3-D)  LiDAR, 
and  more  recently,  flash  3-D  LiDAR  based  on  Geiger  mode  or  linear  mode  avalanche 
photodiode  (i).  For  this  work,  we  focus  on  airborne  scanning  3-D  LiDAR,  which  has  been 
widely  used  for  applications  that  include  sensing  topology  and  bathymetry,  as  well  as  urban 
terrain  mapping. 

Typically,  airborne  scanning  3-D  LiDAR  systems  employ  pulsed  laser  light  that  is  swept 
laterally  by  an  oscillating  scanner  mirror,  sampling  at  a  rate  between  10  and  200  Hz,  depending 
on  the  specific  sensor  system.  As  the  aircraft  moves  along  its  flight  path,  the  linear  scan  produces 
a  Z-shaped  pattern  (2).  Measurements  along  the  scan  lines  are  discretized  to  points,  typically  by 
applying  a  threshold  on  the  intensity,  resulting  in  two  to  six  returns  per  pulse  measurement. 

Points  are  given  a  corrected  geo-position  using  calibrated  latitudinal  and  longitudinal  coordinates 
by  correlating  time-stamped  measurements  from  an  inertial  measurement  unit  and  global 
positioning  system  sensor.  The  measurements  are  based  on  the  time-of-flight  principle,  giving 
each  return  its  distance  from  the  sensor,  and  thus,  allowing  an  elevation  value  to  be  computed. 
The  resulting  scanning  data  constitute  a  non-uniformly  sampled  (with  respect  to  the  x-  and  y- 
spatial  directions)  3-D  point  cloud.  The  points  from  the  Z-shaped  scans  are  not  evenly  spaced, 
but  are  rather  semi-randomly  distributed  due  to  flight  path  perturbations  caused  by  wind,  changes 
in  aircraft  heading  and  velocity,  and  other  mechanical/environmental  factors.  Measured  points 
may  also  overlap  from  adjacent  flight  paths,  resulting  in  data  with  locally  varying  point  densities. 
LiDAR  sensors  are  usually  mounted  in  the  nadir-pointing  orientation,  but  off-nadir  pointing 
angles  may  occur  either  by  design  or  due  to  perturbations  in  the  flight  path.  For  off-nadir 
pointing  angles,  “shadows”  cast  by  tall  objects  or  structures  induce  holes  in  the  data  where  no 
points  were  collected.  Due  to  these  factors,  LiDAR  data  can  exhibit  a  high  degree  of  local 
variability  with  respect  to  the  sample  spacing. 

To  interpret  and  further  process  the  3-D  point  cloud  data,  these  data  are  often  converted  to  digital 
elevation  models  (DEMs).  DEMs  are  created  either  in  raster  format,  where  elevation  information 
is  available  at  every  point  on  a  uniform  grid/raster,  or  as  vector-based  triangulated  irregular 
network  consisting  of  irregularly  distributed  nodes  and  lines  with  3-D  coordinates  arranged  in  a 
network  of  non-overlapping  triangles.  The  focus  of  this  work  is  on  rasterized  DEMs,  which  are 
essential  in  various  applications  such  as  terrain  modeling  and  hydrological  modeling  (3,  4). 
Several  interpolation  techniques  have  been  developed  to  compute  rasterized  DEMs,  including  the 
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triangulated  irregular  network  (TIN)  interpolation  (5),  natural  neighbor  (NN)  (6),  ordinary 
kriging  (7),  and  universal  kriging  (5).  NN  and  TIN  interpolation  techniques  are  computationally 
efficient  and  most  commonly  employed  in  commercial  software.  The  most  prevalent  LiDAR 
exploitation  software  is  the  Quick  Terrain  Modeler  (QT  Modeler)  developed  by  Applied 
Imagery,  which  is  widely  used  by  the  military  and  the  intelligence  community.  To  fill  holes  and 
generate  a  rasterized  DEM  (referred  to  as  a  digital  surface  model  in  the  user’s  manual),  QT 
Modeler  applies  adaptive  triangulation,  which  is  a  modified  TIN  technique.  QT  Modeler’s 
adaptive  triangulation  generates  3-D  triangles  across  empty  regions  and  then  samples  the 
elevation  value  of  the  triangle  surface  at  the  spatial  coordinates  of  the  empty  cell  to  interpolate  its 
elevation  value  (9). 

While  TIN  and  NN  techniques  are  computationally  efficient,  these  techniques  may  result  in 
draping  effects  wherever  shadows  occur  and  near  the  edges  of  objects  like  buildings  and  trees. 
Kriging  methods  are  advanced  geostatistical  algorithms  based  on  regionalized  variable  theory. 
However,  kriging  assumes  that  the  spatial  variation  in  the  z-values  is  statistically  homogeneous 
throughout  the  surface.  While  this  assumption  may  hold  in  natural  scenes  for  applications  like 
soil-landscape  and  hydrological  modeling,  it  is  unlikely  to  be  valid  for  dense  urban  scenes  with 
buildings  of  all  heights,  as  well  as  all  sorts  of  trees  and  vehicles.  Although  LiDAR  data  for  soil- 
landscape  and  hydrological  modeling  continue  to  be  acquired,  LiDAR  data  collections  of  urban 
scenes  have  increasingly  been  gaining  importance  and  attention.  In  order  to  ameliorate  the 
draping  effects  inherent  to  current  interpolation  techniques  for  urban  LiDAR  point  clouds,  we 
need  to  ensure  that  the  edge  information  is  well  preserved  during  the  interpolation  procedure  in 
order  to  compute  the  uniformly  gridded  DEMs. 

In  this  report,  we  propose  a  partial  differential  equations  (PDE)  based  approach  to  interpolate  the 
3-D  point  cloud  onto  a  densely  sampled  (i.e.,  upsampled)  uniform  grid.  First,  all  available 
LiDAR  elevation  data  are  used  to  populate  an  upsampling  grid.  Because  the  sample  spacing  (in 
x-  and  y-  spatial  directions)  in  the  upsampling  grid  is  much  smaller  than  the  average  point  cloud 
sample  spacing,  the  populated  upsampling  grid  is  rather  sparse.  Since  PDEs  have  been  applied  to 
automated  image  interpolation  (10)  in  the  past,  we  postulate  that  the  same  principle  can  be  used 
for  the  interpolation  of  elevation  information  on  a  uniform  grid  from  the  sparsely  populated 
upsampling  grid  by  propagating  the  existing  elevation  data  to  compute  the  rest  of  the  points  on 
the  upsampling  grid.  A  second-order  PDE-like  heat  diffusion  equation  can  diffuse  a  value 
(elevation  information  here)  across  an  area  over  time  (11).  However,  this  approach  would  work 
well  only  when  the  missing  grid  points  are  surrounded  by  homogeneous  regions  (areas  of  same 
elevation).  It  would  fail  when  there  is  a  sharp  edge,  for  instance,  in  the  LiDAR  shadow  regions 
along  the  edges  of  buildings  and  would  lead  to  the  same  draping  effects  observable  in  linear 
interpolation  techniques.  This  problem  is  addressed  in  this  report  by  employing  higher-order 
(e.g.,  fourth-order)  PDEs  (12),  through  which  we  can  maintain  smoothness  in  homogeneous 
regions  while  propagating  the  sharp  edge  information  in  the  scene. 
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The  rest  of  the  report  is  organized  as  follows.  In  section  2,  the  concept  of  PDE-based  image 
inpainting  is  introduced.  Section  3  describes  how  this  image  inpainting  technique  is  applied  to 
perform  uniform  grid  upsampling  of  LiDAR  point  cloud.  Simulation  results  are  presented  in 
section  4.  Finally,  section  5  concludes  the  report  with  some  remarks  about  the  proposed  method. 


2.  PDE  Based  Image  Inpainting 


Traditionally,  image  inpainting  was  done  by  skilled  image  restoration  artists  in  museums  to 
restore  ancient  paintings  that  had  cracks  or  gaps  due  to  aging.  More  recently,  image  inpainting 
has  been  widely  used  in  digital  image  processing  for  applications  such  as  image  restoration  and 
super-resolution.  Image  inpainting  can  be  divided  into  two  main  categories:  exemplar-based 
inpainting  and  PDE-based  inpainting.  Exemplar-based  inpainting  technique  uses  the  whole 
image  to  find  clues  as  to  what  constitutes  the  missing  part  of  the  image  by  checking  for  repetitive 
patterns  and  other  clues.  While  they  are  extremely  useful  to  inpaint  textures  in  the  images,  these 
techniques  require  a  lot  of  training  data  to  learn  what  kind  of  patterns  to  look  for  in  the  images. 
On  the  other  hand,  PDE-based  inpainting  techniques  do  not  use  the  whole  image,  but  just  the 
local  regions  around  the  missing  parts  of  the  image.  Furthermore,  PDE-based  techniques  do  not 
need  any  kind  of  training  process  because  they  do  not  use  any  prior  information  regarding  the 
images.  A  good  review  of  the  PDE  based  inpainting  techniques  was  written  by  Chan  and  Vese 
(13). 

The  main  principle  of  PDE-based  techniques  is  to  minimize  the  energy  in  the  image  with  a 
constraint  that  the  difference  between  the  original  image  and  inpainted  image  in  non-inpainted 
areas  is  within  a  certain  threshold.  For  example,  if  it  denotes  the  original  image  (with  no  holes) 
over  a  2-D  domain  Cl,  it0  denotes  the  observed  image  (with  holes)  over  a  2-D  observed  domain 
D,  D  c  Cl,  then  PDE-based  image  inpainting  algorithms  minimize  the  energy  in  the  image 
denoted  by  E  [it]  over  the  domain  Cl  with  a  constraint  that  it  and  it0  are  as  similar  to  each  other 
as  possible  over  the  domain  D  (13).  Mathematically  it  can  be  written  as, 

minE[it]  such  that  E[u/u0]  <  a2,  (1) 

where  a2  represents  the  white  noise  in  the  observed  image.  This  problem  can  be  written  as  an 
unconstrained  optimization  problem  given  by 

min  E[u]  +  A  E[u/u0\,  (2) 

where  the  Lagrange  multiplier  X  denotes  the  balance  between  energy  minimization  and  image 
fitting.  One  can  observe  from  this  equation  that  PDE-based  image  inpainting  produces  an 
inpainted  image  with  least  energy  and  that  best  fits  the  observed  image.  In  all  PDE-based 
techniques,  the  data  fitting  model  is  considered  to  be 
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(3) 


E[u/u0]  =  j~y  JD  (u  —  u0)2dxdy. 

This  is  equivalent  to  determining  the  inpainted  image  that  best  fits  the  observed  image  in  a  least 
square  error  sense.  The  performance  of  the  inpainting  algorithm  is  crucially  dependent  on  the 
choice  of  the  energy  functional  used  in  equation  2.  In  the  literature,  there  have  been  attempts  to 
define  this  energy  functional  in  different  ways,  which  leads  to  different  heat  flow  PDE  equations. 
The  choice  of  the  energy  functional  essentially  determines  how  the  heat  energy  is  diffused  from 
the  surrounding  regions  into  the  holes  or  missing  parts  of  the  image. 

One  of  the  first  attempts  made  was  to  define  the  energy  functional  E  [it]  using  the  Sobolev  norm 
(14),  i.e.,  E\u]  =  f  |  Vit|2dxdy.  Substituting  this  energy  functional  and  the  data  fitting  model 
defined  in  equation  3  into  equation  2,  one  can  obtain  the  optimization  problem  as 

min  fn  |Vii|2dxdy  +  ( u-u0)2dxdy .  (4) 

The  above  optimization  problem  can  be  solved  by  applying  the  Euler-Lagrange  (E-L)  system  of 
equations  defined  for  the  extremum  (minimum  here)  or  stationary  point  of  energy  functionals, 
such  as  the  one  defined  in  equation  4.  The  definition  and  proof  of  the  existence  of  Euler- 
Lagrange  (E-L)  equations  can  be  found  in  reference  15.  The  E-L  equations  for  the  minimization 
problem  in  equation  4  leads  to  a  Fourier  heat  diffusion  equation,  in  which  the  PDE  is  given  by 

ut  —  A  it  +  A(u  —  it0),  u(t  —  0)  =  it0.  (5) 

Here,  ut  is  the  diffused  image  or  inpainted  image  at  any  time  t,  A  represents  the  Laplacian 
operator,  and  it0  is  the  observed  image.  This  equation  represents  the  diffusion  of  temperature 
across  a  metal  over  time.  In  case  of  image  inpainting,  the  same  equation  may  represent  the 
diffusion  of  gray  values  across  the  image,  i.e.,  from  surrounding  regions  into  the  missing  parts 
and  holes.  However,  this  energy  functional  and  this  PDE  are  effective  only  when  the  holes  are 
present  within  and  surrounded  by  homogeneous  regions.  This  method  also  fails  when  there  are 
sharp  edges  present  around  the  missing  parts  of  the  image. 

In  order  to  deal  with  this  shortcoming,  the  Total  Variational  (TV)  model  was  first  introduced  for 
image  denoising  by  Rudin  et  al.  (11),  and  later,  used  for  image  inpainting  by  Chan  and  Shen 
(10).  In  the  TV  model,  the  energy  functional  is  given  by  the  total  variation  of  u,  that  is  |  D(u)  \ 
over  fl,  and  the  optimization  problem  is 

min  fn  |D(u)|2dxdy  +  ( u-u0)2dxdy .  (6) 

The  E-L  system  of  equations  for  this  problem  yield  an  anisotropic  or  nonlinear  diffusion 
equation  as  shown  in 

Ut  =  V-  ©  +  A(U  ~  U° u &  =  °)  =  uo>  (7) 
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where  V.  represents  the  divergence  operator.  The  TV  model  yields  a  second-order  PDE  and  can 
handle  the  edges  in  the  image  to  a  certain  extent  and  propagate  them  into  the  missing  parts  along 
with  the  homogenous  regions. 

There  are  other  higher-order  PDE  models  in  the  literature,  which  can  handle  edge  information 
along  with  the  diffusion  of  smooth  regions,  including  the  Mumford-Shah  model  (16)  and 
Mumford-Shah-Euler  model  (17).  There  are  two  disadvantages  with  all  these  techniques.  First, 
these  PDE  techniques  cannot  handle  large  gaps  or  holes  in  the  image  if  they  fall  along  the  edges 
and  produce  a  stair-casing  effect  over  smooth  regions.  Second,  they  are  computationally  very 
expensive  to  implement. 

In  order  to  deal  with  the  problem  of  large  holes,  inpainting  of  binary  images  using  modified 
Cahn-Hilliard  PDE  equation  was  proposed  by  Bertozzi  et  al.  (18),  and  later  extended  to  grayscale 
images  by  Burger  et  al.  (12).  This  new  technique  is  called  TV  —  H  -1  inpainting  technique  and 
the  fourth-order  PDE  or  the  diffusion  equation  is  given  by 

ut  =  -A  (v.  -  A(u  -  u0 ),  u(t  =  0)  =  u„.  (8) 

In  order  to  discretize  this  PDE  for  numerical  implementation,  a  convexity  splitting  scheme  is 
used  as  shown  in  reference  12  by  introducing  two  constants,  Cx  and  C2,  where  they  are  set  to 
C1  >  1/e  and  C2  >  7.  The  resulting  discrete  time-stepping  scheme  is  given  by 

+  CtAAuk+1  +  C2uk+1  =  CtAAuk  -  A  (v.  +  C2uk  -  A(uk  -  u0).  (9) 

Here  r  represents  the  time-differential  over  which  the  energy  is  being  diffused.  For  complete 
details  on  the  proof  of  existence  of  solution  and  convergence,  please  refer  to  the  work  of  Burger 
et  al.  (12). 

Based  on  the  value  of  A,  one  can  balance  between  smoothness  of  the  final  image  and  fit  between 
the  final  image  and  the  observed  image.  If  A  is  set  to  a  high  value,  the  fit  between  the  final  image 
and  the  observed  image  is  given  more  importance  and  the  final  image  will  look  very  similar  to 
the  observed  image  in  the  non-missing  regions.  On  the  other  hand,  if  A  is  set  to  a  low  value,  the 
final  image  is  very  smooth  since  the  minimization  of  the  energy  functional  is  given  the  priority. 
This  option  is  useful  when  the  noise  in  the  observed  image  needs  to  be  removed  while  it  is  being 
inpainted  into  the  final  image. 


3.  Application  to  3-D  LiDAR  Point  Cloud 


In  order  to  apply  the  TV  —  H  1  inpainting  technique  to  LiDAR  data,  we  consider  the  (X,  Y ) 
coordinates  in  the  data  as  the  2-D  image  domain,  while  the  elevation  or  Z  coordinate  is 
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considered  equivalent  to  the  gray  levels  in  an  image.  There  are  two  tasks  to  be  performed: 

(1)  upsample  the  LiDAR  point  cloud  and  (2)  rasterize  the  point  cloud.  First,  a  uniform  grid  or 
raster  of  2-D  points  ((X,  Y)  coordinates)  is  defined  over  the  area  of  interest.  The  upsampling  of 
the  data  is  performed  by  choosing  the  distance  between  the  grid  points  to  be  at  least  half  the 
average  distance  between  the  points  in  the  dataset.  Once  the  grid  is  generated,  it  is  populated 
with  the  elevation  information  from  the  original  point  cloud  using  a  nearest  neighbor  selection. 
For  each  datum  in  the  point  cloud,  the  nearest  point  on  the  (X,  Y )  grid  is  determined  in  terms  of 
Euclidean  distance  and  the  elevation  of  the  original  point  is  assigned  to  the  nearest  grid  point  in 
the  raster.  After  assigning  the  elevation  information  of  all  the  cloud  points  to  the  uniform  grid 
points,  many  empty  cells  (i.e.,  holes)  will  be  present  in  the  uniform  grid.  A  real-world  example 
of  a  uniform  grid  after  the  nearest  neighbor  interpolation  is  shown  in  figure  1.  Here  the  white 
points  represent  the  points  of  the  uniform  grid  where  the  elevation  information  is  present.  The 
black  regions  in  the  figure  represent  the  holes.  These  holes  are  filled  by  diffusing  the 
homogeneous  regions  surrounding  them  into  the  holes  and  propagating  edge  information 
wherever  necessary  using  TV  —  H  _1  PDE-based  inpainting  technique. 


Points  in  region  placed  in  upsampling-grid 


Figure  1.  Points  from  3-D  LiDAR  cloud  placed  on  upsampled  uniform  grid. 
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4.  Simulation  Results 


The  airborne  LiDAR  data  used  in  this  work  was  collected  by  Merrick  and  Company  over  the  city 
of  Lubbock,  TX,  using  a  fixed-wing  aircraft.  All  data  and  imagery  shown  in  this  work  are 
unclassified.  In  this  section,  four  sets  of  results  are  presented  to  show  the  effectiveness  of  the 
fourth-order  PDE  based  inpainting  to  upsample  and  rasterize  the  LiDAR  point  cloud.  The 
average  distance  between  the  3-D  points  in  the  original  data  is  0.4  m.  The  12  error  between  the 
uniform  grid  point  cloud  with  holes  and  the  inpainted  uniform  grid  point  cloud  is  used  as  the 
termination  criterion.  When  this  error  falls  below  a  predefined  threshold,  the  algorithm  is 
considered  as  converged  and  the  time-stepping  process  is  then  terminated. 

The  first  dataset  used  for  our  simulation  covers  an  area  around  a  building  (called  building  “A”  in 
this  report),  which  contains  objects  like  buildings,  trees,  and  shrubs.  The  original  point  cloud 
collected  over  the  region  is  shown  in  figure  2a.  The  distance  between  grid  points  is  set  to  0.05  m 
i.e.,  the  data  are  being  upsampled  by  8  times  along  both  x  and  y  axes.  The  balance  parameter  A  is 
set  to  100.  The  upsampled  point  cloud  on  a  uniform  grid  is  illustrated  in  figure  2b.  This  figure 
shows  the  dense  points  on  the  upsampled  grid,  as  well  as  the  elevation  information  obtained  after 
the  inpainting  process.  In  order  to  compare  the  inpainting  technique  to  the  standard  interpolation 
technique  used  in  QT  Modeler,  the  surfaces  generated  by  QT  Modeler  using  standard  internal 
interpolation  of  the  original  point  cloud,  as  well  as  using  the  upsampled  uniform  grid  point 
cloud,  are  presented  in  figures  2c  and  d,  respectively.  One  can  observe  that  the  surface  model 
obtained  using  upsampled  point  cloud  is  sharper  than  the  one  generated  using  original  point 
cloud,  especially  on  the  roof  of  the  buildings,  as  well  as  the  drop  between  the  tree  and  the  shrub 
in  the  region. 
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(b)  Upsampled  inpainted  point  cloud  on  a 
uniform  grid  (using  proposed  technique) 


(d)  Upsampled  inpainted  surface  viewed  in  QT 
Modeler  (using  proposed  technique) 


Figure  2.  Upsampling  and  inpainting  of  building  A  area. 

The  second  and  third  datasets  used  in  this  work  focus  on  another  building  (referred  to  as  building 
“B”).  The  original  point  cloud  over  the  building  B  area  (displayed  in  QT  Modeler  as  a  surface)  is 
shown  in  figure  3.  Two  large  gaps  exist  next  to  the  building:  Gap  1  is  shown  in  figure  3a  and 
Gap  2  in  figure  3b.  These  gaps  might  be  incurred  through  line-of-sight  occlusion  by  the  side  of 
the  building  along  the  aircraft’s  flight  path  or  due  to  the  presence  of  water  that  absorbed  most  of 
the  LiDAR  signal  in  those  areas.  To  illustrate  the  benefit  of  the  upsampling  inpainting  approach, 
sections  of  the  building  that  include  Gap  1  and  Gap  2  were  used  in  our  simulations. 
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Figure  3.  Two  large  gaps  of  building  B. 


The  simulation  results  for  Gap  1  of  building  B  are  presented  in  figure  4.  The  original  point  cloud 
collected  over  Gap  1  is  shown  in  figure  4a.  The  distance  between  grid  points  is  set  to  0.1  m,  i.e., 
the  data  are  being  upsampled  by  4  times  along  both  x  and  y  axes.  The  balance  parameter  A  is  set 
to  100.  The  upsampled  point  cloud  on  a  uniform  grid  is  illustrated  in  figure  4b.  The  surfaces 
generated  by  QT  Modeler  using  standard  internal  interpolation  of  the  original  point  cloud  and 
using  the  upsampled  uniform  grid  point  cloud  are  presented  in  figures  4c  and  d,  respectively. 

One  can  observe  from  figure  4c  that  TIN  interpolation  creates  a  flat  surface  with  constant  slope 
to  connect  the  top  of  the  building  to  the  edge  of  the  gap.  On  the  other  hand,  the  same  wall  in  the 
upsampled  inpainted  surface  in  figure  4d  has  a  much  sharper  interpolated  slope,  especially  near 
the  top  of  the  building.  This  visually  helps  to  improve  the  perception  of  the  wall  of  the  building 
and  proves  that  the  TV  —  H  inpainting  technique  is  able  to  propagate  the  information  from 
the  surrounding  areas  into  the  holes  better  than  simple  interpolation. 
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(b)  Upsampled  inpainted  point  cloud  on  a 
uniform  grid  (using  proposed  technique) 


(c)  Original  surface  viewed  in  QT  Modeler 


A 


(d)  Upsampled  inpainted  surface  viewed  in  QT 
Modeler  (using  proposed  technique) 


Figure  4.  Upsampling  and  inpainting  of  building  B  -  Gap  1. 
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The  simulation  results  for  Gap  2  of  building  B  are  presented  in  figure  5.  The  original  point  cloud 
collected  over  Gap  2  is  shown  in  figure  5a.  The  distance  between  grid  points  is  set  to  0.2  m,  i.e., 
the  data  are  being  upsampled  by  2  times  along  both  x  and  y  axes.  The  balance  parameter  A  is 
again  set  to  100.  The  upsampled  point  cloud  on  a  uniform  grid  is  illustrated  in  figure  5b.  The 
surfaces  generated  by  QT  Modeler  using  standard  internal  interpolation  of  the  original  point 
cloud  and  using  the  upsampled  uniform  grid  point  cloud  are  presented  in  figures  5c  and  5d, 
respectively.  The  piece-wise  linear  interpolation  of  original  point  cloud  produced  spurious  slant 
surfaces  from  the  top  of  the  building  to  the  ground  level  that  cover  the  hole  in  the  image. 
However,  the  PDE-based  inpainting  technique  generates  a  more  vertical  wall  from  the  top  of  the 
building  to  the  ground  level  and  propagates  the  edge  information  from  the  surrounding  shrubs  in 
to  the  hole. 
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(b)  Upsampled  inpainted  point  cloud  on  a 
uniform  grid  (using  proposed  technique) 


(c)  Original  surface  viewed  in  QT  Modeler 


(d)  Upsampled  inpainted  surface  viewed  in  QT 
Modeler  (using  proposed  technique) 


Figure  5.  Upsampling  and  inpainting  of  building  B  -  Gap  2. 
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The  fourth  and  final  dataset  used  in  this  work  covers  the  region  around  a  silo,  which  allows  us  to 
demonstrate  the  effects  of  the  balance  parameter  A.  Using  PDE-based  techniques,  not  only  can 
we  obtain  sharp  edges,  but  we  can  also  smooth  the  noisy  data  points  on  curved  surfaces.  For  this 
experiment,  A  is  set  to  10,  i.e.,  priority  is  being  given  to  the  minimization  of  the  energy 
functional  or  smooth  diffusion  of  the  homogeneous  regions  over  the  fitting  between  the  inpainted 
point  cloud  and  original  point  cloud.  The  results  for  the  silo  region  are  presented  in  figure  6.  The 
original  point  cloud  collected  over  the  silo  region  is  shown  in  figure  6a.  The  distance  between 
grid  points  is  set  to  0.1  m,  i.e.,  the  data  are  being  upsampled  by  4  times  along  both  x  and  y  axes. 
The  upsampled  point  cloud  on  a  uniform  grid  is  illustrated  in  figure  6b.  The  surfaces  generated 
by  QT  Modeler  using  standard  internal  interpolation  of  the  original  point  cloud  and  using  the 
upsampled  uniform  grid  point  cloud  are  presented  in  figures  6c  and  d,  respectively.  It  is  clear  that 
the  TIN  interpolation  of  the  original  point  cloud  generated  jagged  structures  on  the  top  of  the 
silo,  but  the  PDE-based  inpainting  technique  produced  a  smooth  hemispherical  surface,  while 
maintaining  flat  vertical  wall  around  the  silo. 
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(a)  Original  point  cloud 


(b)  Upsampled  inpainted  point  cloud  on  a 
uniform  grid  (using  proposed  technique) 


(c)  Original  surface  viewed  in  QT  Modeler  (d)  Upsampled  inpainted  surface  viewed  in  QT 

Modeler  (using  proposed  technique) 


Figure  6.  Upsampling  and  inpainting  of  a  silo. 


5.  Conclusion 


In  this  report,  a  fourth-order  PDE  based  inpainting  technique  is  proposed  and  evaluated  to 
upsample  3-D  LiDAR  point  cloud  on  to  a  uniform  grid.  This  technique,  called  TV  —  H  -1 
inpainting  technique,  is  based  on  modified  Cahn-Hilliard  gradient  flow  equation.  Our  simulation 
results  have  shown  that  this  technique  reduces  the  draping  effects  that  are  usually  generated  by 
the  standard  piece- wise  linear  interpolation  technique  used  for  upsampling  the  LiDAR  point 
cloud.  This  technique  also  ensures  smoothness  over  homogeneous  regions,  while  propagating 
edges  over  large  holes. 
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Appendix.  Code  for  the  Uniform  Grid  Upsampling  of  3-D  LiDAR  Point  Cloud 
Data 


The  following  is  the  code  used  for  the  uniform  grid  upsampling  of  3-D  LiDAR  point  cloud  data. 

clear  all;close  all;clc 

%  READING  FILE 

%  NOTE:  Easier  to  remove  all  header  information  (.xyz  -  . txt  file),  and  read 

in  numeric  data 

fid=f open (' 504  building.txt'); 

%  read  numeric  data 

points=textscan ( f id,  ' %f  %f  %f  %d  %d  %d  %d  %d  %d  %f  %d  %d  %f'); 

X=points { 1 }  ; 

Y=points { 2  }  ; 

Z=points { 3  }  ; 
clear  points 
f close (fid) ; 

X0=max (X) -min (X)  %  calculating  horizontal  span  of  point  cloud  (m) 

Y0=max ( Y) -min ( Y)  %  calculating  vertical  span  of  point  cloud  (m) 

%  truncating  to  smaller  region  within  point  cloud  using  X  Y  utm  coordinates 
ind=f ind ( (Y>3 . 719447*10 A6) & (Y<3 . 719497*10 A6) & (X>2 . 32715* 10 A5) & (X<2 . 32742*10 A5 
) )  ; 

I=X ( ind) ; 

J=Y ( ind) ; 

E=Z ( ind) ; 

s=5;  %  plotting  subsampled  version  of  point  cloud 

figure 

plot3k ( { X ( 1 : s : end)  Y(l:s:end)  Z ( 1 : s : end) } ) ;  axis  equal; 

set (gca,  ' xtick '  ,  [ ] )  ; 

set (gca, ' ytick '  , [ ] ) ; 

set(gca,  ztick ' , []); 

set (gcf,  ' color ',  [1  1  1]); 

xlabel (  '  X  '  )  ; 

ylabel (  '  Y  '  )  ; 

title ('LiDAR  Point  Cloud'); 
figure 

plot3k({I  J  E});  axis  equal 

set (gca,  'xtick '  ,  [ ] )  ; 

set (gca, 'ytick ' , [ ] ) ; 

set(gca, 'ztick', []); 

set (gcf ,' color ',  [1  1  1]); 

xlabel ( ' X ' ) ; 

ylabel ( ' Y ' ) ; 

title (' Region  of  LiDAR  Point  Cloud'); 


X0=max ( I ) -min ( I ) ; 
Y0=max ( J) -min ( J) ; 
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I vec=min ( I ) : (max ( I ) -min ( I ) ) : max ( I ) ; 

Jvec=min ( J) : (max ( J) -min ( J) ) :max ( J) ; 

%  find  minimum  distance  between  points  (wrt  X  and  Y,  NOT  Z  coordinates) 
%  NOTE:  D  IS  NxN  MATRIX  -  MEMORY  ISSUE  DEPENDING  ON  SIZE  OF  POINT  CLOUD 
%  NICE  TO  FIRST  EXAMINE  DISTRIBUTION  OF  DISTANCES  BETWEEN  ALL  POINTS 
D=L2_distance ( [I  J]  '  ,  [I  J]  '  )  ; 

Dpos=D (D~=0 ) ; 

N=length (Dpos) ; 
minD=min (min (Dpos) ) 

INF=le6 ; 

D(D==0)=INF;  %  set  distance  between  a  point  and  itself  to  be  INF 

minDvec=min (D) ;  %  find  distance  between  every  pair  of  closest  points 
[N  X] =hist (minDvec, 100 ) ; 
figure 
stem (X, N) ; 

xlabel (' Distance  between  every  pair  of  closest  points'); 


S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

oooooooooooooooooooooooooooooooooooo 

%  define  spacing  (minD  or  smaller)  % 

S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

oooooooooooooooooooooooooooooooooooo 

spacing=0.2;  %  spacing  can  be  defined  by  looking  at  histogram  of  using  minD 
LX=max ( I ) -min ( I ) ; 

LY=max ( J) -min ( J) ; 

%  form  grid  for  x  and  y  based  on  minimum  spacing  between  points 
gridXvec=min (I) : spacing: max (I) ; 
gridYvec=min ( J) : spacing :max ( J) ; 

%  making  grid  even 

if  mod (length (gridXvec) , 2) ==1 

gridXvec= [gridXvec  gridXvec (end) fspacing] ; 

end 

if  mod ( length (gridYvec) , 2 ) ==1 

gridYvec= [gridYvec  gridYvec (end) fspacing] ; 

end 

[gridXmatr  gridYmatr] =meshgrid (gridXvec, gridYvec) ; 

[r  c] =size (gridXmatr) ; 

%  let  gridX  and  gridY  be  centers  of  each  cell  in  high  resolution  grid 
gridX=reshape (gridXmatr, 1 , r*c) +spacing/ 2 ; 
gridY=reshape (gridYmatr, 1 , r*c) tspacing/ 2 ; 

D=L2_distance ( [gridX '  gridY ' ]  ' ,  [ I  J]  ' )  ; 

%  find  the  closest  cell  center  to  each  of  the  points 
[val  ind] =min (D,  [  ] , 1 ) ; 
clear  D; 

Zvals= zeros ( 1 , r*c) ; 

Zvals ( ind) =E ; 

I  =  f ind ( Zvals~=0 ) ; 
figure 

plot3k ( { gridX ( I ) '  gridY (I)'  Zvals (I)'});  axis  equal 

set (gca,  ' xtick '  ,  [ ]  )  ; 

set  (gca,  ' ytick ' ,  []); 

set(gca, ' ztick ' , []); 

set (gcf ,  color',  [1  1  1] )  ; 

xlabel ( ' X ' ) ; 

ylabel ( ' Y ' ) ; 
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title (' Region  of  LiDAR  Point  Cloud'); 


X=reshape (gridX, r, c) ; 

Y=reshape (gridY, r, c) ; 

Z=reshape (Zvals, r, c) ; 
mask  =  zeros (size  (Z) ) ; 

I  =  find ( Z==0 ) ; 
ma  s  k ( I )  =  1 ; 
figure 

colormap (gray (256) ) ; 

imagesc ( f lipud ( Z ) ) ;  axis  image; 

set (gca,  ' xtick ' ,  []); 

set (gca,  ytick ' , [ ] ) ; 

title (' Points  in  region  placed  in  upsampling-grid'); 


g=Z  ; 

[m,  n] =size (g) ; 
g=double (g) ; 
gmax  =  max (max (g) ) ; 
g  =  g . /max (max (g) ) ; 
clims= [0  1 ] ; 
figure 

imagesc (g, dims ) ;  axis  image;  axis  off;  colormap (gray) ; 


S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

oooooooooooooooooooo 


Definition  of 


PARAMETERS : 


S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S- 

ooooooooooooooooooooooo 


hl=l; 
h2=l  ; 
dt=100 ; 

epsilon  =  0.01; 
ep2  =  epsilonA2; 
lpower=2 ; 

lambdaO  =  10Alpower; 
cl=l/epsilon; 
c2=lambda0 ; 
approxerr=10 A  (-21)  ; 

%  Definition  of  the  inpainting  domain: 
uO  =  mask; 

lambda  =  ones (m, n) * lambdaO ; 

[j,k]  =  find(u0==l); 
sd  =  size  ( j , 1) ; 
for  i=l : sd 


1 ambda ( j ( i ) , k ( i ) )  =  0; 

end 

clear  j  k  sd 
figure 

imagesc(uO);  axis  image;  axis  off;  colormap (gray) 


%%%%%%  Part  of  this  code  was  adopted  from  the  code  written  by  the  authors 
%%%%%%  of  Ref.  12 

%  Diagonalize  the  Laplace  Operator  by:  Lu  +  uL  =>  D  QuQ  +  QuQ  D, 

%  where  Q  is  nonsingular,  the  matrix  of  eigenvectors  of  L 
%  and  D  is  a  diagonal  matrix.  We  have  to  compute  QuQ. 

%  This  we  can  do  in  a  fast  way  by  using  the  f ft-transform: 
uO  =  g; 

%  Initialization  of  u  and  its  Fourier  transform: 
u=u0  ; 
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ubar=fft2 (u) ; 

Iu0bar=fft2 (lambda . *u)  ; 
lubar=luObar ; 
curv=zeros (m, n) ; 

%Definition  of  the  Laplace  operator  with  periodic  boundary  conditions: 
Ll=l/ (hlA2) * (diag (-2*ones (m, 1) )  +  diag (ones (m-1 , 1 ) , 1 )  +  diag (ones (m-1 , 1 ), - 

1)  .  .  . 

+  diag (ones ( 1 , 1 ), m-1 ) +  diag(ones(l,l),-m+l)); 

L2=l/ (h2A2) * (diag (-2*ones (n, 1) )  +  diag (ones (n-1 , 1 ), 1 )  +  diag (ones (n-1 , 1 ), - 

1)  .  .  . 

+  diag (ones ( 1 , 1 ), n-1 ) +  diag (ones (1, 1) , -n+1) ) ; 

%  Computation  of  the  eigenvalues  of  LI : 
for  j=l :m 

eigvl(j)  =  2*  (cos (2* (j-1) *pi/m) -1) ; 

Lambdal ( j , j ) =eigvl ( j  )  ; 

end 

%  Computation  of  the  eigenvalues  of  L2 : 
for  j  =1 : n 

eigv2(j)  =  2* (cos (2* (j-1) *pi/n) -1) ; 

Lambda2 ( j , j ) =eigv2 ( j )  ; 

end 

Denominator=l /hi A2 . *Lambdal*ones (m, n)  +  l/h2A2 . *ones (m, n) *Lambda2; 

%  Now  we  can  write  the  above  equation  in  much  simpler  way  and  compute 

%  the  solution  ubar 

figure 

it=0  ; 

err=l ; 

tic 

while  err  >  approxerr 


err 

%  Computation  of  the  tv-seminorm: 
%  estimate  derivatives 


ux  = 

(u  ( : ,  [2  :n  n]  ; 

1 -u  (  :  , 

[1  1 : n-1 ] ) 

)  /  2  ; 

uy  = 

(u ( [2 :m  m] , : ; 

l-u([l 

1 :m-l ] , :) 

)  /  2; 

uxx 

=  u  ( : ,  [2  :n  n]  J 

1 +u ( : , 

[1  1 : n-1 ] ) 

-  2  *  u  ; 

uyy 

=  u ( [2 :m  m] , : ; 

I+U([l 

1 :m-l ] , :) 

-  2  *  u  ; 

Dp  = 

u ( [2 :m  m] ,  [2 

:  n  n]) 

+u ( [1  l:m- 

-1] ,  [1 

1 : n-1] )  ; 

Dm  = 

uxy 

u ( [1  1 :m-l ] , 
=  (Dp-Dm) / 4 ; 

[2  :  n  n 

] ) +u  (  [2 :m 

m]  ,  [1 

1 : n-1] )  ; 

%  compute  flow 

Num  =  uxx. * (ep2+uy. A2) -2*ux. *uy. *uxy+uyy. * (ep2+ux. A2) ; 

Den  =  (ep2+ux. A2+uy. A2)  . A  (3/2) ; 
curv  =  Num. /Den; 

lapcurvbar  =  Denominator . *fft2 (curv) ; 

ubar  =  ( (l+c2*dt+dt*cl*Denominator . A2) . *ubar  -  dt*lapcurvbar  ... 

+  dt* ( luObar-lubar ) ) . / (l+dt*c2+dt*cl*Denominator . A2) ; 
unew  =  real (if ft2 (ubar) ) ; 

%error  in  1A2  norm: 

err=  sum ( sum ( (unew-u) . A2 ) ) / (m*n) ; 

u=unew; 

lubar=f f t2 (lambda . *u)  ; 
it=it+l ; 

end 

duration=toc; 

imagesc(u);  axis  image;  axis  off;  colormap (gray) ; 
pause (0.01)  ; 

disp(['u  is  inpainted  in  '  num2str(it)  '  iterations  and  '  num2str (duration) 
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CPU  seconds ' ] ) 


%  Displaying  the  upsampled  3D  point  cloud  on  a  uniform  grid 
u  =  u*gmax/max (max (u) ) ; 

I  =  X; 

J  =  Y; 

E  =  u; 
figure 

plot3k({I  J  E});  axis  equal 
set (gca,  xtick ' ,  [ ]  )  ; 
set  (gca,  ytick '  ,  [ ] )  ; 
set(gca, 'ztick', []); 
set(gcf, ‘color', [1  1  1  ]  )  ; 
xlabel ( ' X ' ) ; 
ylabel ( ' Y ' ) ; 

title (' Region  of  LiDAR  Point  Cloud'); 
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Intentionally  lelt  blank. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


2-D 

two-dimensional 

3-D 

three-dimensional 

DEMs 

digital  elevation  models 

E-L 

Euler-Lagrange 

LiDAR 

light  detection  and  ranging 

NN 

natural  neighbor 

PDEs 

partial  differential  equations 

QT  Modeler 

Quick  Terrain  Modeler 

TIN 

triangulated  irregular  network 

TV 

Total  Variational 

23 


NO.  OF 
COPIES 

1 

(PDF) 


1 

(PDF) 

4 

(PDFS) 


2 

(PDFS) 


2 

(PDFS) 


18 

(PDFS) 


ORGANIZATION 

ADMNSTR 

DEFNS  TECHF  INFO  CTR 
ATTN  DTICOCP 

GOVT  PRINTG  OFC 
A  MAFHOTRA 

US  ARMY  CERDEC  NVESD 
ATTN  L  GRACEFFO 
ATTN  M  GROENERT 
ATTN  J  HILGER 
ATTN  J  WRIGHT 

US  ARMY  AMRDEC 

ATTN  RDMR  WDG  I  J  MILLS 

ATTN  RDMR  WDG  S  D  WAAGEN 

US  ARMY  RSRCH  OFFICE 
ATTN  RDRL  ROI  C  L  DAI 
ATTN  RDRL  ROI  M  J  LAVERY 

US  ARMY  RSRCH  LAB 

ATTN  IMAL  HRA  MAIL  &  RECORDS  MGMT 

ATTN  RDRL  CIO  LL  TECHL  LIB 

ATTN  RDRL  SE  P  PERCONTI 

ATTN  RDRL  SES  J  EICKE 

ATTN  RDRL  SES  M  D’ONOFRIO 

ATTN  RDRL  SES  N  NASRABADI 

ATTN  RDRL  SES  E  R  RAO 

ATTN  RDRL  SES  E  A  CHAN 

ATTN  RDRL  SES  E  H  KWON 

ATTN  RDRL  SES  E  S  YOUNG 

ATTN  RDRL  SES  E  J  DAMMANN 

ATTN  RDRL  SES  E  D  ROSARIO 

ATTN  RDRL  SES  E  H  BRANDT 

ATTN  RDRL  SES  E  S  HU 

ATTN  RDRL  SES  E  M  THIELKE 

ATTN  RDRL  SES  E  P  RAUSS 

ATTN  RDRL  SES  E  P  GURRAM 

ATTN  RDRL  SES  E  C  REALE 


24 


